//
//  CCMatrix4.cpp
//  ShadowMapping
//
//  Created by Daniel Corleone on 16/6/14.
//  Copyright (c) 2014年 iCloverStudio. All rights reserved.
//

#ifndef __CC_MATRIX4_CPP__
#define __CC_MATRIX4_CPP__

#include "CCMatrixs.h"
#include <algorithm>

const Real DEG2RAD = 3.141593 / 180.0;

template<typename T>
Matrix4<T>::Matrix4()
: m00(0), m01(0), m02(0), m03(0)
, m10(0), m11(0), m12(0), m13(0)
, m20(0), m21(0), m22(0), m23(0)
, m30(0), m31(0), m32(0), m33(0)
{}

template<typename T>
Matrix4<T>::Matrix4(value_type _00, value_type _01, value_type _02, value_type _03,
                    value_type _10, value_type _11, value_type _12, value_type _13,
                    value_type _20, value_type _21, value_type _22, value_type _23,
                    value_type _30, value_type _31, value_type _32, value_type _33)
: m00(_00), m01(_01), m02(_02), m03(_03)
, m10(_10), m11(_11), m12(_12), m13(_13)
, m20(_20), m21(_21), m22(_22), m23(_23)
, m30(_30), m31(_31), m32(_32), m33(_33)
{}

template<typename T>
Matrix4<T>::Matrix4(const value_type m[16])
: m00(m[0]), m01(m[1]), m02(m[2]), m03(m[3])
, m10(m[4]), m11(m[5]), m12(m[6]), m13(m[7])
, m20(m[8]), m21(m[9]), m22(m[10]), m23(m[11])
, m30(m[12]), m31(m[13]), m32(m[14]), m33(m[15])
{}

template<typename T>
Matrix4<T>::Matrix4(const Matrix4<T>& other)
: m00(other.m[0]), m01(other.m[1]), m02(other.m[2]), m03(other.m[3])
, m10(other.m[4]), m11(other.m[5]), m12(other.m[6]), m13(other.m[7])
, m20(other.m[8]), m21(other.m[9]), m22(other.m[10]), m23(other.m[11])
, m30(other.m[12]), m31(other.m[13]), m32(other.m[14]), m33(other.m[15])
{}

template<typename T>
CC_INLINE Real
Matrix4<T>::getCofactor(value_type m0, value_type m1, value_type m2, value_type m3,
                        value_type m4, value_type m5, value_type m6, value_type m7,
                        value_type m8, value_type m9, value_type m10, value_type m11,
                        value_type m12, value_type m13, value_type m14, value_type m15)
{
    return (m0 * (m4 * m8 - m5 * m7) -
            m1 * (m3 * m8 - m5 * m6) +
            m2 * (m3 * m7 - m4 * m6));
}

template<typename T>
CC_INLINE void
Matrix4<T>::set(const T m[16])
{
    m00 = m[0]; m01 = m[1]; m02 = m[2]; m03 = m[3];
    m10 = m[4]; m11 = m[5]; m12 = m[6]; m13 = m[7];
    m20 = m[8]; m21 = m[9]; m22 = m[10]; m23 = m[11];
    m30 = m[12]; m31 = m[13]; m32 = m[14]; m33 = m[15];
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::set(value_type _00, value_type _01, value_type _02, value_type _03,
                value_type _10, value_type _11, value_type _12, value_type _13,
                value_type _20, value_type _21, value_type _22, value_type _23,
                value_type _30, value_type _31, value_type _32, value_type _33)
{
    m00 = _00; m01 = _01; m02 = _02; m03 = _03;
    m10 = _10; m11 = _11; m12 = _12; m13 = _13;
    m20 = _20; m21 = _21; m22 = _22; m23 = _23;
    m30 = _30; m31 = _31; m32 = _32; m33 = _33;
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::setRow(Int32 index, const value_type row[4]) {
    Int32 _index = index * 4;
    m[_index] = row[0];
    m[_index + 1] = row[1];
    m[_index + 2] = row[2];
    m[_index + 3] = row[3];
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::setRow(Int32 index, const Vector3<T>& other) {
    Int32 _index = index * 4;
    m[_index] = other.x;
    m[_index + 1] = other.y;
    m[_index + 2] = other.z;
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::setRow(Int32 index, const Vector4<T>& other) {
    Int32 _index = index * 4;
    m[_index] = other.x;
    m[_index + 1] = other.y;
    m[_index + 2] = other.z;
    m[_index + 3] = other.w;
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::setColumn(Int32 index, const value_type col[4]) {
    m[index] = col[0];
    m[index + 4] = col[1];
    m[index + 8] = col[2];
    m[index + 12] = col[3];
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::setColumn(Int32 index, const Vector3<T>& v) {
    m[index] = v.x;
    m[index + 4] = v.y;
    m[index + 8] = v.z;
    return;
}

template<typename T>
CC_INLINE void
Matrix4<T>::setColumn(Int32 index, const Vector4<T>& v) {
    m[index] = v.x;
    m[index + 4] = v.y;
    m[index + 8] = v.z;
    m[index + 12] = v.w;
    return;
}

template<typename T>
CC_INLINE const T*
Matrix4<T>::get() const
{
    return m;
}

template<typename T>
CC_INLINE const T*
Matrix4<T>::getTranspose()
{
    tm[0] = m[0]; tm[1] = m[4]; tm[2] = m[8]; tm[3] = m[12];
    tm[4] = m[1]; tm[5] = m[5]; tm[6] = m[9]; tm[7] = m[13];
    tm[8] = m[2]; tm[9] = m[6]; tm[10] = m[10]; tm[11] = m[14];
    tm[12] = m[3]; tm[13] = m[7]; tm[14] = m[11]; tm[15] = m[15];
    return tm;
}

template<typename T>
CC_INLINE Real
Matrix4<T>::getDeterminant()
{
    return (m[0] * getCofactor(m[5], m[6], m[7], m[9], m[10], m[11], m[13], m[14], m[15]) -
            m[1] * getCofactor(m[4], m[6], m[7], m[8], m[10], m[11], m[12], m[14], m[15]) +
            m[2] * getCofactor(m[4], m[5], m[7], m[8], m[9],  m[11], m[12], m[13], m[15]) -
            m[3] * getCofactor(m[4], m[5], m[6], m[8], m[9],  m[10], m[12], m[13], m[14]));
}


template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::identity()
{
    m[0] = m[5] = m[10] = m[15] = 1.0;
    m[1] = m[2] = m[3] = m[4] = m[6] = m[7] = m[8] = m[9] = m[11] = m[12] = m[13] = m[14] = 0.0;
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::transpose()
{
    std::swap(m[1],  m[4]);
    std::swap(m[2],  m[8]);
    std::swap(m[3],  m[12]);
    std::swap(m[6],  m[9]);
    std::swap(m[7],  m[13]);
    std::swap(m[11], m[14]);
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::invert()
{
    // If the 4th row is [0,0,0,1] then it is affine matrix and
    // it has no projective transformation.
    if(m[12] == 0 && m[13] == 0 && m[14] == 0 && m[15] == 1)
        this->invertAffine();
    else
    {
        this->invertGeneral();
        /*@@ invertProjective() is not optimized (slower than generic one)
         if(fabs(m[0]*m[5] - m[1]*m[4]) > 0.00001f)
         this->invertProjective();   // inverse using matrix partition
         else
         this->invertGeneral();      // generalized inverse
         */
    }
    
    
    return *this;
}

///////////////////////////////////////////////////////////////////////////////
// compute the inverse of 4x4 Euclidean transformation matrix
//
// Euclidean transformation is translation, rotation, and reflection.
// With Euclidean transform, only the position and orientation of the object
// will be changed. Euclidean transform does not change the shape of an object
// (no scaling). Length and angle are reserved.
//
// Use inverseAffine() if the matrix has scale and shear transformation.
//
// M = [ R | T ]
//     [ --+-- ]    (R denotes 3x3 rotation/reflection matrix)
//     [ 0 | 1 ]    (T denotes 1x3 translation matrix)
//
// y = M*x  ->  y = R*x + T  ->  x = R^-1*(y - T)  ->  x = R^T*y - R^T*T
// (R is orthogonal,  R^-1 = R^T)
//
//  [ R | T ]-1    [ R^T | -R^T * T ]    (R denotes 3x3 rotation matrix)
//  [ --+-- ]   =  [ ----+--------- ]    (T denotes 1x3 translation)
//  [ 0 | 1 ]      [  0  |     1    ]    (R^T denotes R-transpose)
///////////////////////////////////////////////////////////////////////////////
template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::invertEuclidean()
{
    // transpose 3x3 rotation matrix part
    // | R^T | 0 |
    // | ----+-- |
    // |  0  | 1 |
    T tmp;
    tmp = m[1];  m[1] = m[4];  m[4] = tmp;
    tmp = m[2];  m[2] = m[8];  m[8] = tmp;
    tmp = m[6];  m[6] = m[9];  m[9] = tmp;
    
    // compute translation part -R^T * T
    // | 0 | -R^T x |
    // | --+------- |
    // | 0 |   0    |
    T x = m[3];
    T y = m[7];
    T z = m[11];
    m[3]  = -(m[0] * x + m[1] * y + m[2] * z);
    m[7]  = -(m[4] * x + m[5] * y + m[6] * z);
    m[11] = -(m[8] * x + m[9] * y + m[10]* z);
    
    // last row should be unchanged (0,0,0,1)
    return *this;
}

///////////////////////////////////////////////////////////////////////////////
// compute the inverse of a 4x4 affine transformation matrix
//
// Affine transformations are generalizations of Euclidean transformations.
// Affine transformation includes translation, rotation, reflection, scaling,
// and shearing. Length and angle are NOT preserved.
// M = [ R | T ]
//     [ --+-- ]    (R denotes 3x3 rotation/scale/shear matrix)
//     [ 0 | 1 ]    (T denotes 1x3 translation matrix)
//
// y = M*x  ->  y = R*x + T  ->  x = R^-1*(y - T)  ->  x = R^-1*y - R^-1*T
//
//  [ R | T ]-1   [ R^-1 | -R^-1 * T ]
//  [ --+-- ]   = [ -----+---------- ]
//  [ 0 | 1 ]     [  0   +     1     ]
///////////////////////////////////////////////////////////////////////////////
template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::invertAffine()
{
    // R^-1
    Matrix3<T> r(m[0],m[1],m[2], m[4],m[5],m[6], m[8],m[9],m[10]);
    r.invert();
    m[0] = r[0];  m[1] = r[1];  m[2] = r[2];
    m[4] = r[3];  m[5] = r[4];  m[6] = r[5];
    m[8] = r[6];  m[9] = r[7];  m[10]= r[8];
    
    // -R^-1 * T
    T x = m[3];
    T y = m[7];
    T z = m[11];
    m[3]  = -(r[0] * x + r[1] * y + r[2] * z);
    m[7]  = -(r[3] * x + r[4] * y + r[5] * z);
    m[11] = -(r[6] * x + r[7] * y + r[8] * z);
    
    // last row should be unchanged (0,0,0,1)
    //m[12] = m[13] = m[14] = 0.0f;
    //m[15] = 1.0f;
    
    return * this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::invertProjective()
{
    // partition
    Matrix2<T> a(m[0], m[1], m[4], m[5]);
    Matrix2<T> b(m[2], m[3], m[6], m[7]);
    Matrix2<T> c(m[8], m[9], m[12], m[13]);
    Matrix2<T> d(m[10], m[11], m[14], m[15]);
    
    // pre-compute repeated parts
    a.invert();             // A^-1
    Matrix2<T> ab = a * b;     // A^-1 * B
    Matrix2<T> ca = c * a;     // C * A^-1
    Matrix2<T> cab = ca * b;   // C * A^-1 * B
    Matrix2<T> dcab = d - cab; // D - C * A^-1 * B
    
    // check determinant if |D - C * A^-1 * B| = 0
    //NOTE: this function assumes det(A) is already checked. if |A|=0 then,
    //      cannot use this function.
    Real determinant = dcab[0] * dcab[3] - dcab[1] * dcab[2];
    if(fabs(determinant) <= FLOAT_MIN) {
        return identity();
    }
    
    // compute D' and -D'
    Matrix2<T> d1 = dcab;      //  (D - C * A^-1 * B)
    d1.invert();            //  (D - C * A^-1 * B)^-1
    Matrix2<T> d2 = -d1;       // -(D - C * A^-1 * B)^-1
    
    // compute C'
    Matrix2<T> c1 = d2 * ca;   // -D' * (C * A^-1)
    
    // compute B'
    Matrix2<T> b1 = ab * d2;   // (A^-1 * B) * -D'
    
    // compute A'
    Matrix2<T> a1 = a - (ab * c1); // A^-1 - (A^-1 * B) * C'
    
    // assemble inverse matrix
    m[0] = a1[0];  m[1] = a1[1];  m[2] = b1[0];  m[3] = b1[1];
    m[4] = a1[2];  m[5] = a1[3];  m[6] = b1[1];  m[2] = b1[3];
    m[8] = c1[0];  m[9] = c1[1];  m[10]= d1[0];  m[11]= d1[1];
    m[12]= c1[2];  m[13]= c1[3];  m[14]= d1[2];  m[15]= d1[3];
    
    return *this;
}


template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::invertGeneral()
{
    // get cofactors of minor matrices
    float cofactor0 = getCofactor(m[5],m[6],m[7], m[9],m[10],m[11], m[13],m[14],m[15]);
    float cofactor1 = getCofactor(m[4],m[6],m[7], m[8],m[10],m[11], m[12],m[14],m[15]);
    float cofactor2 = getCofactor(m[4],m[5],m[7], m[8],m[9], m[11], m[12],m[13],m[15]);
    float cofactor3 = getCofactor(m[4],m[5],m[6], m[8],m[9], m[10], m[12],m[13],m[14]);
    
    // get determinant
    float determinant = m[0] * cofactor0 - m[1] * cofactor1 + m[2] * cofactor2 - m[3] * cofactor3;
    if(::fabs(determinant) <= FLOAT_MIN) {
        return identity();
    }
    
    // get rest of cofactors for adj(M)
    float cofactor4 = getCofactor(m[1],m[2],m[3], m[9],m[10],m[11], m[13],m[14],m[15]);
    float cofactor5 = getCofactor(m[0],m[2],m[3], m[8],m[10],m[11], m[12],m[14],m[15]);
    float cofactor6 = getCofactor(m[0],m[1],m[3], m[8],m[9], m[11], m[12],m[13],m[15]);
    float cofactor7 = getCofactor(m[0],m[1],m[2], m[8],m[9], m[10], m[12],m[13],m[14]);
    
    float cofactor8 = getCofactor(m[1],m[2],m[3], m[5],m[6], m[7],  m[13],m[14],m[15]);
    float cofactor9 = getCofactor(m[0],m[2],m[3], m[4],m[6], m[7],  m[12],m[14],m[15]);
    float cofactor10= getCofactor(m[0],m[1],m[3], m[4],m[5], m[7],  m[12],m[13],m[15]);
    float cofactor11= getCofactor(m[0],m[1],m[2], m[4],m[5], m[6],  m[12],m[13],m[14]);
    
    float cofactor12= getCofactor(m[1],m[2],m[3], m[5],m[6], m[7],  m[9], m[10],m[11]);
    float cofactor13= getCofactor(m[0],m[2],m[3], m[4],m[6], m[7],  m[8], m[10],m[11]);
    float cofactor14= getCofactor(m[0],m[1],m[3], m[4],m[5], m[7],  m[8], m[9], m[11]);
    float cofactor15= getCofactor(m[0],m[1],m[2], m[4],m[5], m[6],  m[8], m[9], m[10]);
    
    // build inverse matrix = adj(M) / det(M)
    // adjugate of M is the transpose of the cofactor matrix of M
    float invDeterminant = 1.0f / determinant;
    m[0] =  invDeterminant * cofactor0;
    m[1] = -invDeterminant * cofactor4;
    m[2] =  invDeterminant * cofactor8;
    m[3] = -invDeterminant * cofactor12;
    
    m[4] = -invDeterminant * cofactor1;
    m[5] =  invDeterminant * cofactor5;
    m[6] = -invDeterminant * cofactor9;
    m[7] =  invDeterminant * cofactor13;
    
    m[8] =  invDeterminant * cofactor2;
    m[9] = -invDeterminant * cofactor6;
    m[10]=  invDeterminant * cofactor10;
    m[11]= -invDeterminant * cofactor14;
    
    m[12]= -invDeterminant * cofactor3;
    m[13]=  invDeterminant * cofactor7;
    m[14]= -invDeterminant * cofactor11;
    m[15]=  invDeterminant * cofactor15;
    
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::translate(const T x, const T y, const T z) {
    m[0] += m[12]*x;   m[1] += m[13]*x;   m[2] += m[14]*x;   m[3] += m[15]*x;
    m[4] += m[12]*y;   m[5] += m[13]*y;   m[6] += m[14]*y;   m[7] += m[15]*y;
    m[8] += m[12]*z;   m[9] += m[13]*z;   m[10]+= m[14]*z;   m[11]+= m[15]*z;
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::translate(const Vector3<T>& other) {
    return translate(other.x, other.y, other.z);
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::rotate(Real angle, const Vector3f& axis) {
    return rotate(angle, axis.x, axis.y, axis.z);
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::rotate(Real angle, const Real x, const Real y, const Real z) {
    Real c = ::cosf(angle * DEG2RAD);    // cosine
    Real s = ::sinf(angle * DEG2RAD);    // sine
    Real xx = x * x;
    Real xy = x * y;
    Real xz = x * z;
    Real yy = y * y;
    Real yz = y * z;
    Real zz = z * z;
    
    // build rotation matrix
    Matrix4<T> m;
    m[0] = xx * (1 - c) + c;
    m[1] = xy * (1 - c) - z * s;
    m[2] = xz * (1 - c) + y * s;
    m[3] = 0;
    m[4] = xy * (1 - c) + z * s;
    m[5] = yy * (1 - c) + c;
    m[6] = yz * (1 - c) - x * s;
    m[7] = 0;
    m[8] = xz * (1 - c) - y * s;
    m[9] = yz * (1 - c) + x * s;
    m[10]= zz * (1 - c) + c;
    m[11]= 0;
    m[12]= 0;
    m[13]= 0;
    m[14]= 0;
    m[15]= 1;
    
    // multiply it
    *this = m * (*this);
    
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::rotateX(Real angle) {
    Real c = ::cosf(angle * DEG2RAD);
    Real s = ::sinf(angle * DEG2RAD);
    Real m4 = m[4], m5 = m[5], m6 = m[6],  m7 = m[7],
    m8 = m[8], m9 = m[9], m10= m[10], m11= m[11];
    
    m[4] = m4 * c + m8 *-s;
    m[5] = m5 * c + m9 *-s;
    m[6] = m6 * c + m10*-s;
    m[7] = m7 * c + m11*-s;
    m[8] = m4 * s + m8 * c;
    m[9] = m5 * s + m9 * c;
    m[10]= m6 * s + m10* c;
    m[11]= m7 * s + m11* c;
    
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::rotateY(Real angle) {
    Real c = ::cosf(angle * DEG2RAD);
    Real s = ::sinf(angle * DEG2RAD);
    Real m0 = m[0], m1 = m[1], m2 = m[2],  m3 = m[3],
    m8 = m[8], m9 = m[9], m10= m[10], m11= m[11];
    
    m[0] = m0 * c + m8 * s;
    m[1] = m1 * c + m9 * s;
    m[2] = m2 * c + m10* s;
    m[3] = m3 * c + m11* s;
    m[8] = m0 *-s + m8 * c;
    m[9] = m1 *-s + m9 * c;
    m[10]= m2 *-s + m10* c;
    m[11]= m3 *-s + m11* c;
    
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::rotateZ(Real angle) {
    Real c = ::cosf(angle * DEG2RAD);
    Real s = ::sinf(angle * DEG2RAD);
    Real m0 = m[0], m1 = m[1], m2 = m[2],  m3 = m[3],
    m4 = m[4], m5 = m[5], m6 = m[6],  m7 = m[7];
    
    m[0] = m0 * c + m4 *-s;
    m[1] = m1 * c + m5 *-s;
    m[2] = m2 * c + m6 *-s;
    m[3] = m3 * c + m7 *-s;
    m[4] = m0 * s + m4 * c;
    m[5] = m1 * s + m5 * c;
    m[6] = m2 * s + m6 * c;
    m[7] = m3 * s + m7 * c;
    
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::scale(const Real s)
{
    return scale(s, s, s);
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::scale(const Real x, const Real y, const Real z)
{
    m[0] = m[0]*x;   m[1] = m[1]*x;   m[2] = m[2]*x;   m[3] = m[3]*x;
    m[4] = m[4]*y;   m[5] = m[5]*y;   m[6] = m[6]*y;   m[7] = m[7]*y;
    m[8] = m[8]*z;   m[9] = m[9]*z;   m[10]= m[10]*z;  m[11]= m[11]*z;
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>
Matrix4<T>::operator+(const Matrix4<T>& other) const {
    return Matrix4<T>(m[0] + other[0], m[1] + other[1], m[2] + other[2], m[3] + other[3],
                      m[4] + other[4], m[5] + other[5], m[6] + other[6], m[7] + other[7],
                      m[8] + other[8], m[9] + other[9], m[10] + other[10], m[11] + other[11],
                      m[12] + other[12], m[13] + other[13], m[14] + other[14], m[15] + other[15]);
}

template<typename T>
CC_INLINE Matrix4<T>
Matrix4<T>::operator-(const Matrix4<T>& other) const {
    return Matrix4<T>(m[0] - other[0], m[1] - other[1], m[2] - other[2], m[3] - other[3],
                      m[4] - other[4], m[5] - other[5], m[6] - other[6], m[7] - other[7],
                      m[8] - other[8], m[9] - other[9], m[10] - other[10], m[11] - other[11],
                      m[12] - other[12], m[13] - other[13], m[14] - other[14], m[15] - other[15]);
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::operator+=(const Matrix4<T>& other) {
    m[0] += other[0]; m[1] += other[1]; m[2] += other[2]; m[3] += other[3];
    m[4] += other[4]; m[5] += other[5]; m[6] += other[6]; m[7] += other[7];
    m[8] += other[8]; m[9] += other[9]; m[10] += other[10]; m[11] += other[11];
    m[12] += other[12]; m[13] += other[13]; m[14] += other[14]; m[15] += other[15];
    return *this;
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::operator-=(const Matrix4<T>& other) {
    m[0] -= other[0]; m[1] -= other[1]; m[2] -= other[2]; m[3] -= other[3];
    m[4] -= other[4]; m[5] -= other[5]; m[6] -= other[6]; m[7] -= other[7];
    m[8] -= other[8]; m[9] -= other[9]; m[10] -= other[10]; m[11] -= other[11];
    m[12] -= other[12]; m[13] -= other[13]; m[14] -= other[14]; m[15] -= other[15];
    return *this;
}

template<typename T>
CC_INLINE Vector4<T>
Matrix4<T>::operator*(const Vector4<T>& other) const {
    return Vector4<T>(m[0]*other.x + m[1]*other.y + m[2]*other.z + m[3]*other.w,
                      m[4]*other.x + m[5]*other.y + m[6]*other.z + m[7]*other.w,
                      m[8]*other.x + m[9]*other.y + m[10]*other.z + m[11]*other.w,
                      m[12]*other.x + m[13]*other.y + m[14]*other.z + m[15]*other.w);
}

template<typename T>
CC_INLINE Vector3<T>
Matrix4<T>::operator*(const Vector3<T>& other) const {
    return Vector3<T>(m[0]*other.x + m[1]*other.y + m[2]*other.z,
                      m[4]*other.x + m[5]*other.y + m[6]*other.z,
                      m[8]*other.x + m[9]*other.y + m[10]*other.z);
}

template<typename T>
CC_INLINE Matrix4<T>
Matrix4<T>::operator*(const Matrix4<T>& other) const {
    return Matrix4<T>();
}

template<typename T>
CC_INLINE Matrix4<T>&
Matrix4<T>::operator*=(const Matrix4<T>& other) {
    *this = *this * other;
    return *this;
}

template<typename T>
CC_INLINE Boolean
Matrix4<T>::operator==(const Matrix4<T>& other) const {
    return ((m[0] == other[0]) && (m[1] == other[1]) && (m[2] == other[2]) && (m[3] == other[3]) &&
            (m[4] == other[4]) && (m[5] == other[5]) && (m[6] == other[6]) && (m[7] == other[7]) &&
            (m[8] == other[8]) && (m[9] == other[9]) && (m[10] == other[10]) && (m[11] == other[11]) &&
            (m[12] == other[12]) && (m[13] == other[13]) && (m[14] == other[14]) && (m[15] == other[15]));
}

template<typename T>
CC_INLINE Boolean
Matrix4<T>::operator!=(const Matrix4<T>& other) const {
    return !(*this == other);
}

template<typename T>
CC_INLINE T
Matrix4<T>::operator[](const Int32 index) const {
    return m[index];
}

template<typename T>
CC_INLINE T&
Matrix4<T>::operator[](const Int32 index) {
    return m[index];
}

template<typename T>
CC_INLINE Matrix4<T>
Matrix4<T>::operator-() const {
    return Matrix4<T>(-m[0], -m[1], -m[2], -m[3],
                      -m[4], -m[5], -m[6], -m[7],
                      -m[8], -m[9], -m[10], -m[11],
                      -m[12], -m[13], -m[14], -m[15]);
}

template<typename T>
CC_INLINE Matrix4<T>
Matrix4<T>::operator*(Real s) const {
    return Matrix4<T>(m[0] * s, m[1] * s, m[2] * s, m[3] * s,
                      m[4] * s, m[5] * s, m[6] * s, m[7] * s,
                      m[8] * s, m[9] * s, m[10] * s, m[11] * s,
                      m[12] * s, m[13] * s, m[14] * s, m[15] * s);
}


#endif  // __CC_MATRIX4_CPP__
